### CREATED BY JESSICA SCHOENEHRR AND NICK WATERBURY ###
### Replication files for "Confessions at the Supreme Court" ###
### File 2 - Replication of Paper Models and Figures ###

### data sources:
### - Confessions of error originally collected by Bruhl (2009)
### 		- Updated by authors through the 2014 term
### - Politicization measure from Patrick Wohlfarth (2009)
###		- Updated by the authors through the 2014 term
### - sgJusticeDistance from the Judicial Common Space (Epstein, Martin, Segal, and Westerland)
###		- http://epstein.wustl.edu/research/JCS.html
### - CSI from Collins and Cooper
###		- https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/UR2KYE
### - SG list of cases and SG as petitioner from Black and Owens (2012)
###		- updated by the authors per the OSG website
###		- https://www.justice.gov/osg/supreme-court-briefs
### - all other data from the Supreme Court DataBase (Spaeth, Epstein, et al.)
### 		- http://scdb.wustl.edu/index.php

library(arm)
library(compactr)
library(stargazer)
library(haven)
library(modelsummary)

library(miceadds)
library(faraway)

library(tidyverse)
library(ggplot2)
library(ggthemes)
library(gridExtra)

####################
####################
####################
### MERITS MODEL ###
####################
####################
####################

data1 <- read_dta("MeritsJusticeCentered20201026.dta")

data1$justice <- as.factor(data1$justice)
data1$votesWithSG <- as.factor(data1$votesWithSG)

#############
### MODEL ###
#############

model4MeritsFullMLM <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase,
							data = data1,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model4MeritsFullMLM$glm_res)
logLik(model4MeritsFullMLM$glm_res)
BIC(model4MeritsFullMLM$glm_res)
nobs(model4MeritsFullMLM$glm_res)

###
# policy change
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4MeritsFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

meritsPolicySeq <- seq(0, 1, length.out = 2)

pMeanMeritsPolicyNoCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pLoMeritsPolicyNoCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pHiMeritsPolicyNoCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))

pMeanMeritsPolicyCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pLoMeritsPolicyCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pHiMeritsPolicyCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))

pMeanMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsPolicySeq))
pLoMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsPolicySeq))
pHiMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsPolicySeq))

for(i in 1:length(meritsPolicySeq)){
	predMeritsPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase
							)
	predMeritsPolicyNoCOE <- as.matrix(predMeritsPolicyNoCOE)
	
	predMeritsPolicyCOE <- cbind(1, 
							1,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase
							)
	predMeritsPolicyCOE <- as.matrix(predMeritsPolicyCOE)
	
	ppMeritsPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppMeritsPolicyCOE <- as.matrix(NA, row = n.draws)
	ppMeritsPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbMeritsPolicyNoCOE <- ilogit(predMeritsPolicyNoCOE %*% t(sim.coef[j, ]))
		ppMeritsPolicyNoCOE[j] <- mean(xbMeritsPolicyNoCOE)
		xbMeritsPolicyCOE <- ilogit(predMeritsPolicyCOE %*% t(sim.coef[j, ]))
		ppMeritsPolicyCOE[j] <- mean(xbMeritsPolicyCOE)
		ppMeritsPolicyCOEDiff[j] <- mean(xbMeritsPolicyCOE - xbMeritsPolicyNoCOE)
	}
	
	pMeanMeritsPolicyNoCOE <- mean(ppMeritsPolicyNoCOE)
	pLoMeritsPolicyNoCOE <- quantile(ppMeritsPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiMeritsPolicyNoCOE <- quantile(ppMeritsPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsPolicyCOE <- mean(ppMeritsPolicyCOE)
	pLoMeritsPolicyCOE <- quantile(ppMeritsPolicyCOE, 0.025, na.rm = TRUE)
	pHiMeritsPolicyCOE <- quantile(ppMeritsPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsPolicyCOEDiff <- mean(ppMeritsPolicyCOEDiff)
	pLoMeritsPolicyCOEDiff <- quantile(ppMeritsPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiMeritsPolicyCOEDiff <- quantile(ppMeritsPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanMeritsPolicyNoCOE
pLoMeritsPolicyNoCOE
pHiMeritsPolicyNoCOE
nameMeritsPolicyNoCOE <- "No Confession of Error"

pMeanMeritsPolicyCOE
pLoMeritsPolicyCOE
pHiMeritsPolicyCOE
nameMeritsPolicyCOE <- "Policy Change"

pMeanMeritsPolicyCOEDiff
pLoMeritsPolicyCOEDiff
pHiMeritsPolicyCOEDiff

meritsBaseline <- data.frame(nameMeritsPolicyNoCOE, pLoMeritsPolicyNoCOE, pMeanMeritsPolicyNoCOE, pHiMeritsPolicyNoCOE)
colnames(meritsBaseline) <- c("div", "lb", "est", "ub")
rownames(meritsBaseline) <- c()

meritsPolicyCOE <- data.frame(nameMeritsPolicyCOE, pLoMeritsPolicyCOE, pMeanMeritsPolicyCOE, pHiMeritsPolicyCOE)
colnames(meritsPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(meritsPolicyCOE) <- c()

###
# noPolicyChange
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4MeritsFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

meritsNoPolicySeq <- seq(0, 1, length.out = 2)

pMeanMeritsNoPolicyNoCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pLoMeritsNoPolicyNoCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pHiMeritsNoPolicyNoCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))

pMeanMeritsNoPolicyCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pLoMeritsNoPolicyCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pHiMeritsNoPolicyCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))

pMeanMeritsNoPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pLoMeritsNoPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pHiMeritsNoPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsNoPolicySeq))

for(i in 1:length(meritsNoPolicySeq)){
	predMeritsNoPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase
							)
	predMeritsNoPolicyNoCOE <- as.matrix(predMeritsNoPolicyNoCOE)
	
	predMeritsNoPolicyCOE <- cbind(1, 
							0,
							1,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase
							)
	predMeritsNoPolicyCOE <- as.matrix(predMeritsNoPolicyCOE)
	
	ppMeritsNoPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppMeritsNoPolicyCOE <- as.matrix(NA, row = n.draws)
	ppMeritsNoPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbMeritsNoPolicyNoCOE <- ilogit(predMeritsNoPolicyNoCOE %*% t(sim.coef[j, ]))
		ppMeritsNoPolicyNoCOE[j] <- mean(xbMeritsNoPolicyNoCOE)
		xbMeritsNoPolicyCOE <- ilogit(predMeritsNoPolicyCOE %*% t(sim.coef[j, ]))
		ppMeritsNoPolicyCOE[j] <- mean(xbMeritsNoPolicyCOE)
		ppMeritsNoPolicyCOEDiff[j] <- mean(xbMeritsNoPolicyCOE - xbMeritsNoPolicyNoCOE)
	}
	
	pMeanMeritsNoPolicyNoCOE <- mean(ppMeritsNoPolicyNoCOE)
	pLoMeritsNoPolicyNoCOE <- quantile(ppMeritsNoPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiMeritsNoPolicyNoCOE <- quantile(ppMeritsNoPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsNoPolicyCOE <- mean(ppMeritsNoPolicyCOE)
	pLoMeritsNoPolicyCOE <- quantile(ppMeritsNoPolicyCOE, 0.025, na.rm = TRUE)
	pHiMeritsNoPolicyCOE <- quantile(ppMeritsNoPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsNoPolicyCOEDiff <- mean(ppMeritsNoPolicyCOEDiff)
	pLoMeritsNoPolicyCOEDiff <- quantile(ppMeritsNoPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiMeritsNoPolicyCOEDiff <- quantile(ppMeritsNoPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanMeritsNoPolicyNoCOE
pLoMeritsNoPolicyNoCOE
pHiMeritsNoPolicyNoCOE
nameMeritsNoPolicyNoCOE <- "No Confession of Error"

pMeanMeritsNoPolicyCOE
pLoMeritsNoPolicyCOE
pHiMeritsNoPolicyCOE
nameMeritsNoPolicyCOE <- "Other Type of Confession"

pMeanMeritsNoPolicyCOEDiff
pLoMeritsNoPolicyCOEDiff
pHiMeritsNoPolicyCOEDiff

meritsNoPolicyCOE <- data.frame(nameMeritsNoPolicyCOE, pLoMeritsNoPolicyCOE, pMeanMeritsNoPolicyCOE, pHiMeritsNoPolicyCOE)
colnames(meritsNoPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(meritsNoPolicyCOE) <- c()

meritsData <- rbind(meritsBaseline, meritsPolicyCOE, meritsNoPolicyCOE)

meritsData$div <- factor(meritsData$div, levels = c("No Confession of Error", "Policy Change", "Other Type of Confession"))

meritsFigure <- (ggplot(meritsData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
  	theme_light() +
 	ggtitle("SG as Party on the Merits") +
  	scale_x_discrete("", labels = c("No\nConfession\nof Error", "Policy\nChange", "Admission\nof\nMistake")) +
  	scale_y_continuous(limits = c(.38, .86), "Probability the Justice Votes with the OSG", breaks = seq(.38, .86, 0.04)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

meritsFigure

###
# coePastTally
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4MeritsFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

meritsPastSeq <- seq(0, 1, length.out = 2)

pMeanMeritsCOE <- as.matrix(NA, nrow = length(meritsPastSeq))
pLoMeritsCOE <- as.matrix(NA, nrow = length(meritsPastSeq))
pHiMeritsCOE <- as.matrix(NA, nrow = length(meritsPastSeq))

pMeanMeritsPastCOE <- as.matrix(NA, nrow = length(meritsPastSeq))
pLoMeritsPastCOE <- as.matrix(NA, nrow = length(meritsPastSeq))
pHiMeritsPastCOE <- as.matrix(NA, nrow = length(meritsPastSeq))

pMeanMeritsPastDiff <- as.matrix(NA, nrow = length(meritsPastSeq))
pLoMeritsPastDiff <- as.matrix(NA, nrow = length(meritsPastSeq))
pHiMeritsPastDiff <- as.matrix(NA, nrow = length(meritsPastSeq))

for(i in 1:length(meritsPastSeq)){
	predMeritsCOE <- cbind(1, 
							1,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase
							)
	predMeritsCOE <- as.matrix(predMeritsCOE)
	
	predMeritsPastCOE <- cbind(1, 
							0,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							1,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase
							)
	predMeritsPastCOE <- as.matrix(predMeritsPastCOE)
	
	ppMeritsCOE <- as.matrix(NA, row = n.draws)
	ppMeritsPastCOE <- as.matrix(NA, row = n.draws)
	ppMeritsPastDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbMeritsCOE <- ilogit(predMeritsCOE %*% t(sim.coef[j, ]))
		ppMeritsCOE[j] <- mean(xbMeritsCOE)
		xbMeritsPastCOE <- ilogit(predMeritsPastCOE %*% t(sim.coef[j, ]))
		ppMeritsPastCOE[j] <- mean(xbMeritsPastCOE)
		ppMeritsPastDiff[j] <- mean(xbMeritsPastCOE - xbMeritsCOE)
	}
	
	pMeanMeritsCOE <- mean(ppMeritsCOE)
	pLoMeritsCOE <- quantile(ppMeritsCOE, 0.025, na.rm = TRUE)
	pHiMeritsCOE <- quantile(ppMeritsCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsPastCOE <- mean(ppMeritsPastCOE)
	pLoMeritsPastCOE <- quantile(ppMeritsPastCOE, 0.025, na.rm = TRUE)
	pHiMeritsPastCOE <- quantile(ppMeritsPastCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsPastDiff <- mean(ppMeritsPastDiff)
	pLoMeritsPastDiff <- quantile(ppMeritsPastDiff, 0.025, na.rm = TRUE)
	pHiMeritsPastDiff <- quantile(ppMeritsPastDiff, 0.975, na.rm = TRUE)
}


pMeanMeritsCOE
pLoMeritsCOE
pHiMeritsCOE
nameMeritsCOE <- "Policy Change in Current Session"

pMeanMeritsPastCOE
pLoMeritsPastCOE
pHiMeritsPastCOE
nameMeritsPastCOE <- "Confession of Error in Last Session"

pMeanMeritsPastDiff
pLoMeritsPastDiff
pHiMeritsPastDiff

meritsCOE <- data.frame(nameMeritsCOE, pLoMeritsCOE, pMeanMeritsCOE, pHiMeritsCOE)
colnames(meritsCOE) <- c("div", "lb", "est", "ub")
rownames(meritsCOE) <- c()

meritsPastCOE <- data.frame(nameMeritsPastCOE, pLoMeritsPastCOE, pMeanMeritsPastCOE, pHiMeritsPastCOE)
colnames(meritsPastCOE) <- c("div", "lb", "est", "ub")
rownames(meritsPastCOE) <- c()

meritsPastData <- rbind(meritsCOE, meritsPastCOE)

meritsPastData$div <- factor(meritsPastData$div, levels = c("Policy Change in Current Session", "Confession of Error in Last Session"))

meritsPastFigure <- (ggplot(meritsPastData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
	geom_hline(aes(yintercept = 0.64), linetype = "dashed") +
  	theme_light() +
 	ggtitle("SG as Party on the Merits") +
  	scale_x_discrete("", labels = c("Policy Change\nConfession\nin Current Session", "Confession of Error\nin Last Session")) +
  	scale_y_continuous(limits = c(.38, .90), "Probability the Justice Votes with the OSG", breaks = seq(.38, .90, 0.05)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

meritsPastFigure

###################
###################
###################
### CVSG MODELS ###
###################
###################
###################

data2 <- read_dta("CVSGJusticeCentered20201026.dta")

data2$justice <- as.factor(data2$justice)
data2$votesWithSG <- as.factor(data2$votesWithSG)

#############
### MODEL ###
#############

model4CVSGFullMLM <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase,
							data = data2,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model4CVSGFullMLM$glm_res)
logLik(model4CVSGFullMLM$glm_res)
BIC(model4CVSGFullMLM$glm_res)
nobs(model4CVSGFullMLM$glm_res)

###
# policy change
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4CVSGFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

CVSGPolicySeq <- seq(0, 1, length.out = 2)

pMeanCVSGPolicyNoCOE <- as.matrix(NA, nrow = length(CVSGPolicySeq))
pLoCVSGPolicyNoCOE <- as.matrix(NA, nrow = length(CVSGPolicySeq))
pHiCVSGPolicyNoCOE <- as.matrix(NA, nrow = length(CVSGPolicySeq))

pMeanCVSGPolicyCOE <- as.matrix(NA, nrow = length(CVSGPolicySeq))
pLoCVSGPolicyCOE <- as.matrix(NA, nrow = length(CVSGPolicySeq))
pHiCVSGPolicyCOE <- as.matrix(NA, nrow = length(CVSGPolicySeq))

pMeanCVSGPolicyCOEDiff <- as.matrix(NA, nrow = length(CVSGPolicySeq))
pLoCVSGPolicyCOEDiff <- as.matrix(NA, nrow = length(CVSGPolicySeq))
pHiCVSGPolicyCOEDiff <- as.matrix(NA, nrow = length(CVSGPolicySeq))

for(i in 1:length(CVSGPolicySeq)){
	predCVSGPolicyNoCOE <- cbind(1, 
							0,
							0,
							data2$coeCriminalCase,
							data2$coePastTally,
							data2$coeInPrevious,
							data2$politicization,
							data2$sgJusticeDist,
							data2$sgPetitioner,
							data2$constitutionalCase,
							data2$declareUncon,
							data2$CSI,
							data2$realOutlierSignal,
							data2$civilRightsCase
							)
	predCVSGPolicyNoCOE <- as.matrix(predCVSGPolicyNoCOE)
	
	predCVSGPolicyCOE <- cbind(1, 
							1,
							0,
							data2$coeCriminalCase,
							data2$coePastTally,
							data2$coeInPrevious,
							data2$politicization,
							data2$sgJusticeDist,
							data2$sgPetitioner,
							data2$constitutionalCase,
							data2$declareUncon,
							data2$CSI,
							data2$realOutlierSignal,
							data2$civilRightsCase
							)
	predCVSGPolicyCOE <- as.matrix(predCVSGPolicyCOE)
	
	ppCVSGPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppCVSGPolicyCOE <- as.matrix(NA, row = n.draws)
	ppCVSGPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbCVSGPolicyNoCOE <- ilogit(predCVSGPolicyNoCOE %*% t(sim.coef[j, ]))
		ppCVSGPolicyNoCOE[j] <- mean(xbCVSGPolicyNoCOE)
		xbCVSGPolicyCOE <- ilogit(predCVSGPolicyCOE %*% t(sim.coef[j, ]))
		ppCVSGPolicyCOE[j] <- mean(xbCVSGPolicyCOE)
		ppCVSGPolicyCOEDiff[j] <- mean(xbCVSGPolicyCOE - xbCVSGPolicyNoCOE)
	}
	
	pMeanCVSGPolicyNoCOE <- mean(ppCVSGPolicyNoCOE)
	pLoCVSGPolicyNoCOE <- quantile(ppCVSGPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiCVSGPolicyNoCOE <- quantile(ppCVSGPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanCVSGPolicyCOE <- mean(ppCVSGPolicyCOE)
	pLoCVSGPolicyCOE <- quantile(ppCVSGPolicyCOE, 0.025, na.rm = TRUE)
	pHiCVSGPolicyCOE <- quantile(ppCVSGPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanCVSGPolicyCOEDiff <- mean(ppCVSGPolicyCOEDiff)
	pLoCVSGPolicyCOEDiff <- quantile(ppCVSGPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiCVSGPolicyCOEDiff <- quantile(ppCVSGPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanCVSGPolicyNoCOE
pLoCVSGPolicyNoCOE
pHiCVSGPolicyNoCOE
nameCVSGPolicyNoCOE <- "No Confession of Error"

pMeanCVSGPolicyCOE
pLoCVSGPolicyCOE
pHiCVSGPolicyCOE
nameCVSGPolicyCOE <- "Policy Change"

pMeanCVSGPolicyCOEDiff
pLoCVSGPolicyCOEDiff
pHiCVSGPolicyCOEDiff

CVSGBaseline <- data.frame(nameCVSGPolicyNoCOE, pLoCVSGPolicyNoCOE, pMeanCVSGPolicyNoCOE, pHiCVSGPolicyNoCOE)
colnames(CVSGBaseline) <- c("div", "lb", "est", "ub")
rownames(CVSGBaseline) <- c()

CVSGPolicyCOE <- data.frame(nameCVSGPolicyCOE, pLoCVSGPolicyCOE, pMeanCVSGPolicyCOE, pHiCVSGPolicyCOE)
colnames(CVSGPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(CVSGPolicyCOE) <- c()


###
# noPolicyChange
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4CVSGFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

CVSGNoPolicySeq <- seq(0, 1, length.out = 2)

pMeanCVSGNoPolicyNoCOE <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))
pLoCVSGNoPolicyNoCOE <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))
pHiCVSGNoPolicyNoCOE <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))

pMeanCVSGNoPolicyCOE <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))
pLoCVSGNoPolicyCOE <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))
pHiCVSGNoPolicyCOE <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))

pMeanCVSGNoPolicyCOEDiff <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))
pLoCVSGNoPolicyCOEDiff <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))
pHiCVSGNoPolicyCOEDiff <- as.matrix(NA, nrow = length(CVSGNoPolicySeq))

for(i in 1:length(CVSGNoPolicySeq)){
	predCVSGNoPolicyNoCOE <- cbind(1, 
							0,
							0,
							data2$coeCriminalCase,
							data2$coePastTally,
							data2$coeInPrevious,
							data2$politicization,
							data2$sgJusticeDist,
							data2$sgPetitioner,
							data2$constitutionalCase,
							data2$declareUncon,
							data2$CSI,
							data2$realOutlierSignal,
							data2$civilRightsCase
							)
	predCVSGNoPolicyNoCOE <- as.matrix(predCVSGNoPolicyNoCOE)
	
	predCVSGNoPolicyCOE <- cbind(1, 
							0,
							1,
							data2$coeCriminalCase,
							data2$coePastTally,
							data2$coeInPrevious,
							data2$politicization,
							data2$sgJusticeDist,
							data2$sgPetitioner,
							data2$constitutionalCase,
							data2$declareUncon,
							data2$CSI,
							data2$realOutlierSignal,
							data2$civilRightsCase
							)
	predCVSGNoPolicyCOE <- as.matrix(predCVSGNoPolicyCOE)
	
	ppCVSGNoPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppCVSGNoPolicyCOE <- as.matrix(NA, row = n.draws)
	ppCVSGNoPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbCVSGNoPolicyNoCOE <- ilogit(predCVSGNoPolicyNoCOE %*% t(sim.coef[j, ]))
		ppCVSGNoPolicyNoCOE[j] <- mean(xbCVSGNoPolicyNoCOE)
		xbCVSGNoPolicyCOE <- ilogit(predCVSGNoPolicyCOE %*% t(sim.coef[j, ]))
		ppCVSGNoPolicyCOE[j] <- mean(xbCVSGNoPolicyCOE)
		ppCVSGNoPolicyCOEDiff[j] <- mean(xbCVSGNoPolicyCOE - xbCVSGNoPolicyNoCOE)
	}
	
	pMeanCVSGNoPolicyNoCOE <- mean(ppCVSGNoPolicyNoCOE)
	pLoCVSGNoPolicyNoCOE <- quantile(ppCVSGNoPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiCVSGNoPolicyNoCOE <- quantile(ppCVSGNoPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanCVSGNoPolicyCOE <- mean(ppCVSGNoPolicyCOE)
	pLoCVSGNoPolicyCOE <- quantile(ppCVSGNoPolicyCOE, 0.025, na.rm = TRUE)
	pHiCVSGNoPolicyCOE <- quantile(ppCVSGNoPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanCVSGNoPolicyCOEDiff <- mean(ppCVSGNoPolicyCOEDiff)
	pLoCVSGNoPolicyCOEDiff <- quantile(ppCVSGNoPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiCVSGNoPolicyCOEDiff <- quantile(ppCVSGNoPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanCVSGNoPolicyNoCOE
pLoCVSGNoPolicyNoCOE
pHiCVSGNoPolicyNoCOE
nameCVSGNoPolicyNoCOE <- "No Confession of Error"

pMeanCVSGNoPolicyCOE
pLoCVSGNoPolicyCOE
pHiCVSGNoPolicyCOE
nameCVSGNoPolicyCOE <- "Other Type of Confession"

pMeanCVSGNoPolicyCOEDiff
pLoCVSGNoPolicyCOEDiff
pHiCVSGNoPolicyCOEDiff

CVSGNoPolicyCOE <- data.frame(nameCVSGNoPolicyCOE, pLoCVSGNoPolicyCOE, pMeanCVSGNoPolicyCOE, pHiCVSGNoPolicyCOE)
colnames(CVSGNoPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(CVSGNoPolicyCOE) <- c()

CVSGData <- rbind(CVSGBaseline, CVSGPolicyCOE, CVSGNoPolicyCOE)

CVSGData$div <- factor(CVSGData$div, levels = c("No Confession of Error", "Policy Change", "Other Type of Confession"))

CVSGFigure <- (ggplot(CVSGData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
  	theme_light() +
 	ggtitle("SG Responding to CVSG") +
  	scale_x_discrete("", labels = c("No\nConfession\nof Error", "Policy\nChange", "Admission\nof\nMistake")) +
  	scale_y_continuous(limits = c(.38, .86), "Probability the Justice Sides with the OSG", breaks = seq(.38, .86, 0.04)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

CVSGFigure

###
# coePastTally
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4CVSGFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

CVSGPastSeq <- seq(0, 1, length.out = 2)

pMeanCVSGCOE <- as.matrix(NA, nrow = length(CVSGPastSeq))
pLoCVSGCOE <- as.matrix(NA, nrow = length(CVSGPastSeq))
pHiCVSGCOE <- as.matrix(NA, nrow = length(CVSGPastSeq))

pMeanCVSGPastCOE <- as.matrix(NA, nrow = length(CVSGPastSeq))
pLoCVSGPastCOE <- as.matrix(NA, nrow = length(CVSGPastSeq))
pHiCVSGPastCOE <- as.matrix(NA, nrow = length(CVSGPastSeq))

pMeanCVSGPastDiff <- as.matrix(NA, nrow = length(CVSGPastSeq))
pLoCVSGPastDiff <- as.matrix(NA, nrow = length(CVSGPastSeq))
pHiCVSGPastDiff <- as.matrix(NA, nrow = length(CVSGPastSeq))

for(i in 1:length(CVSGPastSeq)){
	predCVSGCOE <- cbind(1, 
							1,
							0,
							data2$coeCriminalCase,
							data2$coePastTally,
							data2$coeInPrevious,
							data2$politicization,
							data2$sgJusticeDist,
							data2$sgPetitioner,
							data2$constitutionalCase,
							data2$declareUncon,
							data2$CSI,
							data2$realOutlierSignal,
							data2$civilRightsCase
							)
	predCVSGCOE <- as.matrix(predCVSGCOE)
	
	predCVSGPastCOE <- cbind(1, 
							0,
							0,
							data2$coeCriminalCase,
							data2$coePastTally,
							1,
							data2$politicization,
							data2$sgJusticeDist,
							data2$sgPetitioner,
							data2$constitutionalCase,
							data2$declareUncon,
							data2$CSI,
							data2$realOutlierSignal,
							data2$civilRightsCase
							)
	predCVSGPastCOE <- as.matrix(predCVSGPastCOE)
	
	ppCVSGCOE <- as.matrix(NA, row = n.draws)
	ppCVSGPastCOE <- as.matrix(NA, row = n.draws)
	ppCVSGPastDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbCVSGCOE <- ilogit(predCVSGCOE %*% t(sim.coef[j, ]))
		ppCVSGCOE[j] <- mean(xbCVSGCOE)
		xbCVSGPastCOE <- ilogit(predCVSGPastCOE %*% t(sim.coef[j, ]))
		ppCVSGPastCOE[j] <- mean(xbCVSGPastCOE)
		ppCVSGPastDiff[j] <- mean(xbCVSGPastCOE - xbCVSGCOE)
	}
	
	pMeanCVSGCOE <- mean(ppCVSGCOE)
	pLoCVSGCOE <- quantile(ppCVSGCOE, 0.025, na.rm = TRUE)
	pHiCVSGCOE <- quantile(ppCVSGCOE, 0.975, na.rm = TRUE)
	
	pMeanCVSGPastCOE <- mean(ppCVSGPastCOE)
	pLoCVSGPastCOE <- quantile(ppCVSGPastCOE, 0.025, na.rm = TRUE)
	pHiCVSGPastCOE <- quantile(ppCVSGPastCOE, 0.975, na.rm = TRUE)
	
	pMeanCVSGPastDiff <- mean(ppCVSGPastDiff)
	pLoCVSGPastDiff <- quantile(ppCVSGPastDiff, 0.025, na.rm = TRUE)
	pHiCVSGPastDiff <- quantile(ppCVSGPastDiff, 0.975, na.rm = TRUE)
}


pMeanCVSGCOE
pLoCVSGCOE
pHiCVSGCOE
nameCVSGCOE <- "Policy Change in Current Session"

pMeanCVSGPastCOE
pLoCVSGPastCOE
pHiCVSGPastCOE
nameCVSGPastCOE <- "Confession of Error in Last Session"

pMeanCVSGPastDiff
pLoCVSGPastDiff
pHiCVSGPastDiff

CVSGCOE <- data.frame(nameCVSGCOE, pLoCVSGCOE, pMeanCVSGCOE, pHiCVSGCOE)
colnames(CVSGCOE) <- c("div", "lb", "est", "ub")
rownames(CVSGCOE) <- c()

CVSGPastCOE <- data.frame(nameCVSGPastCOE, pLoCVSGPastCOE, pMeanCVSGPastCOE, pHiCVSGPastCOE)
colnames(CVSGPastCOE) <- c("div", "lb", "est", "ub")
rownames(CVSGPastCOE) <- c()

CVSGPastData <- rbind(CVSGCOE, CVSGPastCOE)

CVSGPastData$div <- factor(CVSGPastData$div, levels = c("Policy Change in Current Session", "Confession of Error in Last Session"))

CVSGPastFigure <- (ggplot(CVSGPastData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
	geom_hline(aes(yintercept = 0.75), linetype = "dashed") +
  	theme_light() +
 	ggtitle("SG Responding to CVSG") +
  	scale_x_discrete("", labels = c("Policy Change\nConfession\nin Current Session", "Confession of Error\nin Last Session")) +
  	scale_y_continuous(limits = c(.38, .90), "Probability the Justice Sides with the OSG", breaks = seq(.38, .90, 0.05)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

CVSGPastFigure

##################
##################
##################
### VACB MODEL ###
##################
##################
##################

data3 <- read_dta("VoluntaryACBJusticeCentered20201026.dta")

data3$justice <- as.factor(data3$justice)
data3$votesWithSG <- as.factor(data3$votesWithSG)

#############
### MODEL ###
#############

model4VACBFullMLM <- glm.cluster(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase,
							data = data3,
							cluster = "justice",
							family = binomial(link = logit))
							
summary(model4VACBFullMLM$glm_res)
logLik(model4VACBFullMLM$glm_res)
BIC(model4VACBFullMLM$glm_res)
nobs(model4VACBFullMLM$glm_res)

###
# policy change
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4VACBFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

VACBPolicySeq <- seq(0, 1, length.out = 2)

pMeanVACBPolicyNoCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pLoVACBPolicyNoCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pHiVACBPolicyNoCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))

pMeanVACBPolicyCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pLoVACBPolicyCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pHiVACBPolicyCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))

pMeanMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBPolicySeq))
pLoMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBPolicySeq))
pHiMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBPolicySeq))

for(i in 1:length(VACBPolicySeq)){
	predVACBPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase
							)
	predVACBPolicyNoCOE <- as.matrix(predVACBPolicyNoCOE)
	
	predVACBPolicyCOE <- cbind(1, 
							1,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase
							)
	predVACBPolicyCOE <- as.matrix(predVACBPolicyCOE)
	
	ppVACBPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppVACBPolicyCOE <- as.matrix(NA, row = n.draws)
	ppVACBPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbVACBPolicyNoCOE <- ilogit(predVACBPolicyNoCOE %*% t(sim.coef[j, ]))
		ppVACBPolicyNoCOE[j] <- mean(xbVACBPolicyNoCOE)
		xbVACBPolicyCOE <- ilogit(predVACBPolicyCOE %*% t(sim.coef[j, ]))
		ppVACBPolicyCOE[j] <- mean(xbVACBPolicyCOE)
		ppVACBPolicyCOEDiff[j] <- mean(xbVACBPolicyCOE - xbVACBPolicyNoCOE)
	}
	
	pMeanVACBPolicyNoCOE <- mean(ppVACBPolicyNoCOE)
	pLoVACBPolicyNoCOE <- quantile(ppVACBPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiVACBPolicyNoCOE <- quantile(ppVACBPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBPolicyCOE <- mean(ppVACBPolicyCOE)
	pLoVACBPolicyCOE <- quantile(ppVACBPolicyCOE, 0.025, na.rm = TRUE)
	pHiVACBPolicyCOE <- quantile(ppVACBPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBPolicyCOEDiff <- mean(ppVACBPolicyCOEDiff)
	pLoVACBPolicyCOEDiff <- quantile(ppVACBPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiVACBPolicyCOEDiff <- quantile(ppVACBPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanVACBPolicyNoCOE
pLoVACBPolicyNoCOE
pHiVACBPolicyNoCOE
nameVACBPolicyNoCOE <- "No Confession of Error"

pMeanVACBPolicyCOE
pLoVACBPolicyCOE
pHiVACBPolicyCOE
nameVACBPolicyCOE <- "Policy Change"

pMeanVACBPolicyCOEDiff
pLoVACBPolicyCOEDiff
pHiVACBPolicyCOEDiff

VACBBaseline <- data.frame(nameVACBPolicyNoCOE, pLoVACBPolicyNoCOE, pMeanVACBPolicyNoCOE, pHiVACBPolicyNoCOE)
colnames(VACBBaseline) <- c("div", "lb", "est", "ub")
rownames(VACBBaseline) <- c()

VACBPolicyCOE <- data.frame(nameVACBPolicyCOE, pLoVACBPolicyCOE, pMeanVACBPolicyCOE, pHiVACBPolicyCOE)
colnames(VACBPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(VACBPolicyCOE) <- c()

###
# noPolicyChange
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4VACBFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

VACBNoPolicySeq <- seq(0, 1, length.out = 2)

pMeanVACBNoPolicyNoCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pLoVACBNoPolicyNoCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pHiVACBNoPolicyNoCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))

pMeanVACBNoPolicyCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pLoVACBNoPolicyCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pHiVACBNoPolicyCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))

pMeanVACBNoPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pLoVACBNoPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pHiVACBNoPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBNoPolicySeq))

for(i in 1:length(VACBNoPolicySeq)){
	predVACBNoPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase
							)
	predVACBNoPolicyNoCOE <- as.matrix(predVACBNoPolicyNoCOE)
	
	predVACBNoPolicyCOE <- cbind(1, 
							0,
							1,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase
							)
	predVACBNoPolicyCOE <- as.matrix(predVACBNoPolicyCOE)
	
	ppVACBNoPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppVACBNoPolicyCOE <- as.matrix(NA, row = n.draws)
	ppVACBNoPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbVACBNoPolicyNoCOE <- ilogit(predVACBNoPolicyNoCOE %*% t(sim.coef[j, ]))
		ppVACBNoPolicyNoCOE[j] <- mean(xbVACBNoPolicyNoCOE)
		xbVACBNoPolicyCOE <- ilogit(predVACBNoPolicyCOE %*% t(sim.coef[j, ]))
		ppVACBNoPolicyCOE[j] <- mean(xbVACBNoPolicyCOE)
		ppVACBNoPolicyCOEDiff[j] <- mean(xbVACBNoPolicyCOE - xbVACBNoPolicyNoCOE)
	}
	
	pMeanVACBNoPolicyNoCOE <- mean(ppVACBNoPolicyNoCOE)
	pLoVACBNoPolicyNoCOE <- quantile(ppVACBNoPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiVACBNoPolicyNoCOE <- quantile(ppVACBNoPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBNoPolicyCOE <- mean(ppVACBNoPolicyCOE)
	pLoVACBNoPolicyCOE <- quantile(ppVACBNoPolicyCOE, 0.025, na.rm = TRUE)
	pHiVACBNoPolicyCOE <- quantile(ppVACBNoPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBNoPolicyCOEDiff <- mean(ppVACBNoPolicyCOEDiff)
	pLoVACBNoPolicyCOEDiff <- quantile(ppVACBNoPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiVACBNoPolicyCOEDiff <- quantile(ppVACBNoPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanVACBNoPolicyNoCOE
pLoVACBNoPolicyNoCOE
pHiVACBNoPolicyNoCOE
nameVACBNoPolicyNoCOE <- "No Confession of Error"

pMeanVACBNoPolicyCOE
pLoVACBNoPolicyCOE
pHiVACBNoPolicyCOE
nameVACBNoPolicyCOE <- "Other Type of Confession"

pMeanVACBNoPolicyCOEDiff
pLoVACBNoPolicyCOEDiff
pHiVACBNoPolicyCOEDiff

VACBNoPolicyCOE <- data.frame(nameVACBNoPolicyCOE, pLoVACBNoPolicyCOE, pMeanVACBNoPolicyCOE, pHiVACBNoPolicyCOE)
colnames(VACBNoPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(VACBNoPolicyCOE) <- c()

VACBData <- rbind(VACBBaseline, VACBPolicyCOE, VACBNoPolicyCOE)

VACBData$div <- factor(VACBData$div, levels = c("No Confession of Error", "Policy Change", "Other Type of Confession"))

VACBFigure <- (ggplot(VACBData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
  	theme_light() +
 	ggtitle("SG as Voluntary Amicus") +
  	scale_x_discrete("", labels = c("No\nConfession\nof Error", "Policy\nChange", "Admission\nof\nMistake")) +
  	scale_y_continuous(limits = c(.38, .86), "Probability the Justice Sides with the OSG", breaks = seq(.38, .86, 0.04)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

VACBFigure

###
# coePastTally
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4VACBFullMLM$glm_res, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

VACBPastSeq <- seq(0, 1, length.out = 2)

pMeanVACBCOE <- as.matrix(NA, nrow = length(VACBPastSeq))
pLoVACBCOE <- as.matrix(NA, nrow = length(VACBPastSeq))
pHiVACBCOE <- as.matrix(NA, nrow = length(VACBPastSeq))

pMeanVACBPastCOE <- as.matrix(NA, nrow = length(VACBPastSeq))
pLoVACBPastCOE <- as.matrix(NA, nrow = length(VACBPastSeq))
pHiVACBPastCOE <- as.matrix(NA, nrow = length(VACBPastSeq))

pMeanVACBPastDiff <- as.matrix(NA, nrow = length(VACBPastSeq))
pLoVACBPastDiff <- as.matrix(NA, nrow = length(VACBPastSeq))
pHiVACBPastDiff <- as.matrix(NA, nrow = length(VACBPastSeq))

for(i in 1:length(VACBPastSeq)){
	predVACBCOE <- cbind(1, 
							1,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase
							)
	predVACBCOE <- as.matrix(predVACBCOE)
	
	predVACBPastCOE <- cbind(1, 
							0,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							1,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase
							)
	predVACBPastCOE <- as.matrix(predVACBPastCOE)
	
	ppVACBCOE <- as.matrix(NA, row = n.draws)
	ppVACBPastCOE <- as.matrix(NA, row = n.draws)
	ppVACBPastDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbVACBCOE <- ilogit(predVACBCOE %*% t(sim.coef[j, ]))
		ppVACBCOE[j] <- mean(xbVACBCOE)
		xbVACBPastCOE <- ilogit(predVACBPastCOE %*% t(sim.coef[j, ]))
		ppVACBPastCOE[j] <- mean(xbVACBPastCOE)
		ppVACBPastDiff[j] <- mean(xbVACBPastCOE - xbVACBCOE)
	}
	
	pMeanVACBCOE <- mean(ppVACBCOE)
	pLoVACBCOE <- quantile(ppVACBCOE, 0.025, na.rm = TRUE)
	pHiVACBCOE <- quantile(ppVACBCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBPastCOE <- mean(ppVACBPastCOE)
	pLoVACBPastCOE <- quantile(ppVACBPastCOE, 0.025, na.rm = TRUE)
	pHiVACBPastCOE <- quantile(ppVACBPastCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBPastDiff <- mean(ppVACBPastDiff)
	pLoVACBPastDiff <- quantile(ppVACBPastDiff, 0.025, na.rm = TRUE)
	pHiVACBPastDiff <- quantile(ppVACBPastDiff, 0.975, na.rm = TRUE)
}


pMeanVACBCOE
pLoVACBCOE
pHiVACBCOE
nameVACBCOE <- "Policy Change in Current Session"

pMeanVACBPastCOE
pLoVACBPastCOE
pHiVACBPastCOE
nameVACBPastCOE <- "Confession of Error in Last Session"

pMeanVACBPastDiff
pLoVACBPastDiff
pHiVACBPastDiff

VACBCOE <- data.frame(nameVACBCOE, pLoVACBCOE, pMeanVACBCOE, pHiVACBCOE)
colnames(VACBCOE) <- c("div", "lb", "est", "ub")
rownames(VACBCOE) <- c()

VACBPastCOE <- data.frame(nameVACBPastCOE, pLoVACBPastCOE, pMeanVACBPastCOE, pHiVACBPastCOE)
colnames(VACBPastCOE) <- c("div", "lb", "est", "ub")
rownames(VACBPastCOE) <- c()

VACBPastData <- rbind(VACBCOE, VACBPastCOE)

VACBPastData$div <- factor(VACBPastData$div, levels = c("Policy Change in Current Session", "Confession of Error in Last Session"))

VACBPastFigure <- (ggplot(VACBPastData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
	geom_hline(aes(yintercept = 0.68), linetype = "dashed") +
  	theme_light() +
 	ggtitle("SG as Voluntary Amicus") +
  	scale_x_discrete("", labels = c("Policy Change\nConfession\nin Current Session", "Confession of Error\nin Last Session")) +
  	scale_y_continuous(limits = c(.38, .90), "Probability the Justice Sides with the OSG", breaks = seq(.38, .90, 0.05)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

VACBPastFigure

### FIGURES ###

fullBaseline <- grid.arrange(meritsFigure, CVSGFigure, VACBFigure, ncol = 3)

fullPast <- grid.arrange(meritsPastFigure, CVSGPastFigure, VACBPastFigure, ncol = 3)






















